#front matter
rm(list=ls())
library(spdep)
library(mgcv)
library(arm)
library(R2OpenBUGS)
set.seed(314159265)
options(scipen=12)

#set working directory
#setwd('C:/Temp')

# Set-up areal structure using North Carolina example in spdep library 
example(nc.sids) 
n <- length(ncCR85.nb) 
N<-c(n)
coords <- coordinates(nc.sids) 
cards <- card(ncCR85.nb) 
col.w <- nb2listw(ncCR85.nb, style="W") #one weight matrix
col.w.alt <- nb2listw(ncCR85.nb, style="B") #another weight matrix
D<-array(unlist(lapply(col.w.alt$weights,sum))) #D matrix

#global parameters
ridge<- 0.97
constant<- 0
beta1<-.7
beta2<-.5

#local parameter and number of treatments
psi<-seq(.1,.9,.1)
M<-length(psi)

#number of replicate samples
S<-150

#number of MCMC iterations
mcmc<-100000
burn<-10000

#output files
psi.car<-matrix(NA,nrow=S,ncol=M)
beta1.car<-matrix(NA,nrow=S,ncol=M)
beta2.car<-matrix(NA,nrow=S,ncol=M)
beta1.lm<-matrix(NA,nrow=S,ncol=M)
beta2.lm<-matrix(NA,nrow=S,ncol=M)
beta1.car.mae<-matrix(NA,nrow=S,ncol=M)
beta2.car.mae<-matrix(NA,nrow=S,ncol=M)
beta1.lm.mae<-matrix(NA,nrow=S,ncol=M)
beta2.lm.mae<-matrix(NA,nrow=S,ncol=M)
     
#initial values
Initials.Lin<-function(){list(beta=c(.1,.1),constant=.1)}
Initials.CAR<-function(){list(beta=c(.1,.1),constant=c(.1),
phi=c(0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0))}

#GeoBUGS Data#
#num<-c(3,3,5,1,3,3,3,5,4,4,5,5,5,4,3,7,3,8,4,3,2,5,5,7,7,7,6,5,5,5,7,3,5,6,4,6,7,3,9,5,4,7,7,5,2,6,6,8,6,6,7,6,7,7,6,1,6,5,4,3,6,6,7,5,6,4,9,5,5,7,7,5,3,6,4,3,2,5,7,2,4,6,5,4,4,5,4,6,6,2,6,4,4,5,3,5,7,4,2,3)
#adj<-c(2,18,19,1,3,18,2,10,18,23,25,7,6,16,28,5,8,28,4,8,17,6,7,17,20,21,15,16,24,31,3,12,25,26,12,14,26,27,29,10,11,25,26,27,14,15,24,30,37,11,13,29,30,9,13,24,5,9,24,28,31,33,36,7,8,20,1,2,3,19,23,34,39,41,1,18,22,34,8,17,21,8,20,19,32,34,43,46,3,18,25,39,40,9,13,15,16,31,37,54,3,10,12,23,26,40,42,10,11,12,25,27,42,47,11,12,26,29,47,48,5,6,16,36,44,11,14,27,30,48,13,14,29,37,48,9,16,24,33,37,49,54,22,35,46,16,31,36,49,51,18,19,22,41,43,52,32,38,46,53,16,28,33,44,51,57,13,24,30,31,48,54,63,35,53,55,18,23,40,41,50,52,65,68,69,23,25,39,42,50,18,34,39,52,25,26,40,47,50,70,71,22,34,46,52,61,64,65,28,36,45,57,87,44,87,22,32,35,43,53,61,26,27,42,48,67,70,27,29,30,37,47,60,63,67,31,33,51,54,59,62,39,40,42,69,70,71,33,36,49,57,59,74,91,34,39,41,43,64,65,35,38,46,55,61,72,75,24,31,37,49,62,63,79,38,53,58,66,72,75,87,36,44,51,80,87,91,55,66,73,78,81,49,51,62,74,48,63,67,43,46,53,64,72,77,49,54,59,74,79,88,37,48,54,60,67,79,82,43,52,61,65,76,39,43,52,64,68,76,55,58,75,78,47,48,60,63,70,82,86,89,92,39,65,69,76,84,39,50,68,71,84,42,47,50,67,71,85,89,42,50,69,70,84,85,89,53,55,61,75,77,58,78,81,51,59,62,83,88,91,53,55,66,72,64,65,68,61,72,58,66,73,81,90,54,62,63,82,88,96,97,57,91,58,73,78,90,63,67,79,86,94,96,74,88,91,93,95,68,69,71,85,70,71,84,89,67,82,89,92,94,44,45,56,57,62,74,79,83,93,97,67,70,71,85,86,92,78,81,51,57,74,80,83,95,67,86,89,94,83,88,95,97,82,86,92,96,98,83,91,93,79,82,94,97,98,79,88,93,96,98,99,100,94,96,97,100,97,100,97,98,99)
#weights<-c(1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1) 

no.co<-length(ncCR85.nb)

num<-rep(NA,no.co)
for(i in 1:no.co){num[i]<-length(ncCR85.nb[[i]])}

adj<-c(NA)
for(i in 1:no.co){adj<-append(adj,ncCR85.nb[[i]])}
adj<-adj[-1]

weights<-rep(1,length(adj))

##########################  Model Specification  ##########################
#set up loop
M<-length(psi)
j<-1
while(j<=M){
	for(i in 1:S){
	#simulate variables
	x1<-rnorm(n,sd=1)
	x2<-rnorm(n,sd=1)
	tot<-psi[j]^2+(1-psi[j])^2
	phi.0 <- rnorm(n, sd=1.5526751*psi[j]/sqrt(tot)) #clustering standard deviation defined as a multiple of heterogeneous standard deviation (which itself is 1)
	theta <- rnorm(n, sd=(1-psi[j])/sqrt(tot))#does this work better?

	#create spatial error and dependent variable
	DrW <- diag(D) - ridge*listw2mat(col.w.alt) # I - rho * W
	CARcov<-solve(DrW)
	CARcovU = chol(CARcov) #  from  (D-rhoW)^ -1 
	phi<-t(CARcovU)%*%phi.0
	y.star<-constant+beta1*x1+beta2*x2+phi+theta
	y<-ifelse(y.star>0,1,0)

	## PIVOT TO R2OpenBUGS ##
	mydata.car<-bugs.data(list(y=as.numeric(y),x1=as.numeric(x1),x2=as.numeric(x2),N=n,adj=as.numeric(adj),weights=as.numeric(weights),num=as.numeric(num)))
	mod.indep<-bugs(data=mydata.car,inits=Initials.Lin,parameters.to.save=c("beta"),model.file="logitMCindep.txt", n.chains=1, n.burnin=burn, n.iter=mcmc)
	mod.car<-bugs(data=mydata.car,inits=Initials.CAR,parameters.to.save=c("beta","psi"),model.file="logitMCcar.txt", n.chains=1, n.burnin=burn, n.iter=mcmc)
	
	##SAVE RESULTS##
	#Independent Model
	beta1.lm[i,j]<-mod.indep$median$beta[1]
	beta2.lm[i,j]<-mod.indep$median$beta[2]

	#CAR Model
	beta1.car[i,j]<-mod.car$median$beta[1]
	beta2.car[i,j]<-mod.car$median$beta[2]
	psi.car[i,j]<-mod.car$median$psi

	#record absolute errors
	beta1.car.mae[i,j]<-abs(beta1.car[i,j]-beta1)
	beta2.car.mae[i,j]<-abs(beta2.car[i,j]-beta2)
	beta1.lm.mae[i,j]<-abs(beta1.lm[i,j]-beta1)
	beta2.lm.mae[i,j]<-abs(beta2.lm[i,j]-beta2)
	}
	
	#Cycle through outer loop
	cat("Round", j, "Complete \n")
	j<-j+1
}

#create summaries
beta1.lm.sum<-apply(beta1.lm,2,mean)
beta1.car.sum<-apply(beta1.car,2,mean)
beta2.lm.sum<-apply(beta2.lm,2,mean)
beta2.car.sum<-apply(beta2.car,2,mean)
psi.car.sum<-apply(psi.car,2,mean)
mae.beta1.lm.sum<-apply(beta1.lm.mae,2,mean)
mae.beta1.car.sum<-apply(beta1.car.mae,2,mean)
mae.beta2.lm.sum<-apply(beta2.lm.mae,2,mean)
mae.beta2.car.sum<-apply(beta2.car.mae,2,mean)

pdf("psiEst.pdf")
plot(x=psi,y=psi.car.sum)
dev.off()

#Graph the results 
pdf("beta1logitMAEind.pdf")
par(mar=c(5.1,4.5,4.1,1.9))
plot(x=psi,y=mae.beta1.lm.sum,type='h',xlab=expression(psi),ylab=expression(paste("Mean Absolute Error for ", beta[1])),col='blue',lty=2,lwd=2,xlim=c(0,1),ylim=c(0,1.5),axes=F)
axis(1,at=seq(0,1,.1));axis(2);box()
abline(h=0,col='gray60')
par(new=T)
plot(x=psi+.01, y=mae.beta1.car.sum, type='h',xlab="",ylab="",axes=F, col='red', lty=1,lwd=2,xlim=c(0,1),ylim=c(0,1.5))
legend(x=0,y=1.5,legend=c("Simple Linear","CAR"), lty=c(2,1),lwd=c(2,2),col=c('blue','red'))
dev.off()

pdf("beta2logitMAEind.pdf")
par(mar=c(5.1,4.5,4.1,1.9))
plot(x=psi,y=mae.beta2.lm.sum,type='h',xlab=expression(psi),ylab=expression(paste("Mean Absolute Error for ", beta[2])),col='blue',lty=2,lwd=2,xlim=c(0,1),ylim=c(0,1.5),axes=F)
axis(1,at=seq(0,1,.1));axis(2);box()
abline(h=0,col='gray60')
par(new=T)
plot(x=psi+.01, y=mae.beta2.car.sum, type='h',xlab="",ylab="",axes=F, col='red', lty=1,lwd=2,xlim=c(0,1),ylim=c(0,1.5))
legend(x=0,y=1.5,legend=c("Simple Linear","CAR"), lty=c(2,1),lwd=c(2,2),col=c('blue','red'))
dev.off()

pdf("beta1logitInd.pdf")
par(mar=c(5.1,4.5,4.1,1.9))
plot(x=psi,y=beta1.lm.sum,type='l',xlab=expression(psi),ylab=expression(hat(beta[1])),col='blue',lty=2,lwd=2,ylim=c(0,2))
lines(x=psi, y=beta1.car.sum, col='red', lty=3,lwd=2)
abline(h=beta1,col='black',lwd=3)
legend(x=.1,y=.85*beta1,legend=c("Simple Linear","CAR","Population"), lty=c(2,3,1),lwd=c(2,2,3),col=c('blue','red','black'))
dev.off()

pdf("beta2logitInd.pdf")
par(mar=c(5.1,4.5,4.1,1.9))
plot(x=psi,y=beta2.lm.sum,type='l',xlab=expression(psi),ylab=expression(hat(beta[2])),col='blue',lty=2,lwd=2,ylim=c(0,2))
lines(x=psi, y=beta2.car.sum, col='red', lty=3,lwd=2)
abline(h=beta2,col='black',lwd=3)
legend(x=.1,y=.85*beta2,legend=c("Simple Linear","CAR","Population"), lty=c(2,3,1),lwd=c(2,2,3),col=c('blue','red','black'))
dev.off()
